Completion of IGN data, road example

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(happign) # retrieve pnr boundary
#> Warning in CPL_gdal_init(): GDAL Message 1: Unable to find driver DODS to
#> unload from GDAL_SKIP environment variable.

#> Warning in CPL_gdal_init(): GDAL Message 1: Unable to find driver DODS to
#> unload from GDAL_SKIP environment variable.
#> Please make sure you have an internet connection.
#> Use happign::get_last_news() to display latest geoservice news.
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(tmap);tmap_mode("view") # interactive map
#> tmap mode set to interactive viewing
library(osmextract)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
#> Check the package website, https://docs.ropensci.org/osmextract/, for more details.
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

Context

To compare the completeness of IGN data, OpenStreetMap data will be used. They will be retrieved with the package osmextract.

The study area is the territory regional natural park of the Massif des Bauges located in the Northern Alps, in the Bauges massif. It has an area of about 900 km2, spread over the departments of Savoie and Haute-Savoie. It is made up of 67 communes.

Extract PNR boundary

To extract IGN data with happign, it is necessary to have a geographical input. I chose the commune of Annecy located in the PNR of the Massif des Bauges. happign contains the cog_2022 dataset which allows to find the insee code of this commune to retrieve boundary with get_apicarto_commune.

code_insee <- cog_2022[cog_2022$LIBELLE == "Annecy", 1][1] #use cog_2022 happign dataset to find insee code
annecy <- get_apicarto_commune(code_insee) # get annecy boundary

pnr <- get_wfs(annecy, 
               apikey = "environnement",
               layer_name = "PROTECTEDAREAS.PNR:pnr") # search for pnr containinng annecy
qtm(pnr) #quick check

Extract roads data

## From ign bdcarto with happign
all_road_ign <- get_wfs(pnr,
                        apikey = "cartovecto",
                        layer_name = "BDCARTO_BDD_WLD_WGS84G:troncon_route")

road_ign <- st_intersection(all_road_ign, pnr)

## from osm data with osmextract
all_road_osm <- oe_get("Savoie",
                       provider = "openstreetmap_fr",
                       stringsAsFactors = FALSE,
                       force_download = TRUE, 
                       force_vectortranslate = TRUE,
                       quiet = TRUE)

road_osm <- st_intersection(all_road_osm, pnr) |>
   filter(highway %in% c("primary", "primary_link", "secondary","secondary_link", 
                         "tertiary", "tertiary_link", "trunk", "trunk_link", 
                         "residential", "cycleway", "living_street", "unclassified", 
                         "motorway", "motorway_link", "pedestrian", "steps",
                         "track","service", NA))

Results

The data set provided by osmextract contains 4 times more linear than the BDCARTO of the ign. Be careful however, the osmextract data requires more verifications than the IGN data : 25% of the OSM dataset are unqualified roads (highway == NA).

# who has the longest ?
sum(st_length(road_ign))
#> 1683201 [m]
sum(st_length(road_osm))
#> 6773427 [m]

# plot results
map <- tm_shape(pnr)+
   tm_borders()+
tm_shape(road_osm)+
   tm_lines(col = "firebrick", lwd = 2)+
tm_shape(road_ign)+
   tm_lines(col = "#95C019", lwd = 2)+
tm_add_legend(type = "fill",
              labels = c("IGN Roads", "OSM Roads"),
              col = c("#95C019", "firebrick"))

map